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RANDOM  FIELD  SOLUTIONS  INCLUDING  BOUNDARY  CONDITION  UNCERTAINTY 
EOR  THE  STEADY-STATE  GENERALIZED  BURGERS  EQUATION 

LUC  HUYSE*  AND  ROBERT  W.  WALTERSt 

Abstract.  CFD  results  are  subject  to  considerable  uncertainty  associated  with  the  operating  conditions.  Even 
when  the  operational  uncertainty  is  omitted  under  very  controlled  circumstances  during  wind  tunnel  experiments, 
substantial  disagreement  between  experimental  and  CFD  results  persists.  This  discrepancy  must  be  attributed  to  model 
uncertainty.  This  report  discusses  the  various  sources  of  uncertainty.  The  need  for  advanced  uncertainty  modeling  is 
illustrated  by  means  of  a  computationally  inexpensive  1-D  Burgers  equation  model.  We  specifically  address  the 
uncertainty  due  to  missing  variables  (inexact  or  incomplete  differential  equations).  To  this  extent  a  random  field 
model  is  used  for  the  viscosity  and  the  fundamental  differences  between  the  solutions  of  the  stochastic  differential 
equations  and  a  simple  random  variable  model  is  highlighted.  The  Burgers  equation  theoretically  needs  to  be  integrated 
over  an  infinite  domain.  In  a  deterministic  approach,  the  integration  domain  is  cut  off  at  some  far  field  boundary. 
This  truncation  effectively  ignores  all  variability  outside  this  far  field  boundary.  We  present  a  practical  treatment  for 
the  uncertainty  on  the  boundary  conditions.  The  results  indicate  that  ignoring  the  boundary  condition  uncertainty 
dramatically  underestimates  the  variance  of  the  velocity  u{x)  in  the  interior  of  the  domain. 

Key  words,  model  uncertainty,  random  field,  uncertain  boundary  conditions 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction  &  Motivation.  Both  computational  results  and  experimental  measurements  in  aerodynamic 
applications  are  subject  to  considerable  uncertainty  as  illustrated  by  the  scatter  in  the  solutions  submitted  by  various 
researchers  for  a  recent  AIAA  drag  prediction  workshop  [1].  For  these  reasons,  the  current  certification  process  relies 
heavily  on  full-scale  flight  testing.  Dramatic  savings  would  be  achieved  if  the  need  for  expensive  full-scale  flight  tests 
could  be  reduced  through  improved  reliability  (or  dependability)  of  CFD  results. 

The  interest  in  uncertainty  modeling  within  the  CFD  community  has  grown  considerably  in  recent  years.  In  a 
deterministic  approach,  the  performance  of  a  design  is  typically  assessed  for  a  limited  number  of  design  or  operating 
conditions.  Uncertainty  associated  with  the  operating  conditions  (variable  payload,  atmospheric  conditions)  and  fluc¬ 
tuations  in  the  geometry  and  smoothness  of  the  skin  (manufacturing  uncertainty)  may  have  substantial  impact  on  the 
results.  Techniques  for  propagating  these  uncertainties  require  additional  computational  effort  but  are  well  established, 
an  example  application  is  described  in  [30].  Exact  and  approximate  uncertainty  assessment  methods  are  applied  to  an 
airfoil  optimization  problem  in  [17],  where  a  dramatic  improvement  of  the  robustness  of  the  design  was  achieved. 

However,  even  if  these  outside  sources  of  “operational”  uncertainty  are  omitted  under  very  controlled  circum¬ 
stances  during  wind  tunnel  experiments,  substantial  disagreement  between  experimental  and  CFD  results  persists  [1]. 
Assuming  that  the  CFD  model  and  the  experiment  simulate  the  same  thing,  this  discrepancy  must  be  attributed  to 
model  uncertainty.  In  this  paper  we  will  address  the  issue  of  modeling  uncertainty  for  a  computationally  inexpensive 
1-D  Burgers  equation  model. 
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2.  Outline  of  Paper.  After  an  initial  discussion  of  the  different  sources  of  uncertainty  and  quantitative  uncertainty 
modeling  techniques  in  Sections  3  and  4,  the  Burgers  equation  is  introduced  in  Section  5.  Results  for  a  random 
variable  model  of  the  viscosity  are  presented  in  Section  6.  A  random  variable  model  of  the  viscosity  assumes  that 
the  viscosity  remains  constant  throughout  the  entire  domain.  However,  the  temperature  dependence  of  the  viscosity 
introduces  additional  uncertainty  in  the  viscosity  which  is  best  captured  using  a  random  held  model.  This  model 
uncertainty  effectively  transforms  the  Burgers  equation  into  a  stochastic  differential  equation.  Second-order  random 
held  models  are  discussed  in  Sections  7.  Two  discretization  techniques,  midpoint  and  locally  averaged  discretizations, 
are  presented  in  Sections  8  and  9.  Results  for  the  random  held  model  are  given  in  Section  10.  The  boundary  conditions 
for  the  Burgers  equation  are  specihed  at  inhnity.  In  numerical  solution  methods,  the  integration  domain  is  limited  and 
cut  off  at  the  “far  held”  boundary.  The  additional  uncertainty  due  to  this  boundary  condition  treatment  is  assessed  in 
Section  11.  Final  conclusions  are  drawn  in  Section  12. 

3.  Sources  of  Uncertainty.  At  this  point  in  the  discussion  it  is  appropriate  to  spend  some  time  on  nomenclature, 
especially  the  difference  between  “error”,  “scatter”  and  “uncertainty”.  Unfortunately,  these  terms  have  been  used  to 
denote  somewhat  different  things  in  different  helds  of  application.  Historically,  the  CFD  community  has  sometimes 
used  the  term  “uncertainty”  to  denote  “error”  in  relation  to  grid  convergence  studies  [32].  This  confounds  the  issue. 
We  adopt  the  following  dehnitions,  which  are  commonly  used  by  statisticians: 

•  Error  is  a  deterministic  concept  and  is  dehned  as  the  difference  between  the  true  or  exact  answer  to  a  problem 
and  the  answer,  computed  or  measured  using  a  faulty  or  simplihed  theory. 

•  Scatter  measures  the  range  or  spread  of  the  data,  but  gives  no  information  about  the  potential  bias  due  to 
systematic  measurement  error  or  due  to  missing  terms  in  the  CFD  code. 

•  Uncertainty  indicates  that  the  result  can  be  only  known  with  a  limited  amount  of  confidence  for  a  given  level 
of  precision.  This  uncertainty  is  an  inherent  property  of  the  measurement  technique  or  model  description  and 
is  due  to  lack  of  knowledge. 

The  distinction  between  the  deterministic  and  stochastic  nature  is  apparent  in  the  AIAA  dehnitions  as  well  [2]: 

•  Error  is  a  recognizable  dehciency  in  any  phase  or  activity  of  modeling  and  simulation  that  is  not  due  to  lack 
of  knowledge. 

•  Uncertainty  is  a  potential  dehciency  in  any  phase  or  activity  of  modeling  and  simulation  that  is  due  to  lack  of 
knowledge. 

Since  error  is  a  recognizable  dehciency,  all  errors  are  -  in  principle  at  least  -  correctable  and  therefore  deterministic. 
Since  uncertainty  is  caused  by  a  fundamental  lack  of  knowledge,  it  cannot  be  eliminated.  If  a  higher  conhdence  level 
(or  level  of  credibility)  of  the  prediction  is  required,  the  result  can  only  be  given  with  less  precision.  Uncertainty  forces 
a  fundamental  trade-off  between  conhdence  and  precision. 

Two  sources  of  uncertainty  can  be  distinguished  [6], [28]:  inherent  (aleatoric)  and  model  (or  epistemic)  uncer¬ 
tainty. 

•  Inherent  Uncertainty:  this  applies  to  phenomena  which  we  accept  to  be  intrinsically  variable  in  our  model.  A 
descriptive,  often  statistical,  formulation  is  used  for  this  uncertainty  and  no  attempt  is  made  to  explain  any  of 
this  variability.  In  a  mathematical  model,  these  uncertainties  are  typically  described  using  (joint)  probability 
density  functions.  Techniques  for  the  efficient  propagation  of  these  distributions  through  a  model  will  be 
briehy  discussed  in  this  report. 

•  Model  Uncertainty:  three  types  can  be  distinguished  here.  The  hrst  two  types  can  be  caused  by  either  “errors 
of  ignorance”  or  “errors  of  simplihcation”.  Note  that  use  of  the  word  “error”  means  that  these  uncertainties 
originate  in  either  an  acknowledged  deliberate  simplihcation  or  an  acknowledged  lack  of  understanding. 
However,  since  the  exact  solution  is  unknown  is  this  case,  the  errors  can  not  be  corrected  and  we  have  no 
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alternative  but  to  treat  them  in  a  non-deterministic  manner  as  uncertainties. 

1 .  missing  variables:  a  typical  example  is  the  exclusion  of  coupling  effects  from  the  analysis.  For  instance 
the  viscosity  depends  on  the  temperature,  but  we  may  decide  not  to  do  a  full  coupled  aero-thermal 
analysis  and  ignore  the  effect  of  the  temperature  on  the  viscosity.  Ignoring  structural  deformations  of 
the  wing  in  an  aerodynamic  analysis  is  another  example  thereof. 

2.  inexact  functional  form:  quite  often  we  use  linear  equations  whereas  the  original  physics  model  contains 
non-linear  terms.  The  linearization  introduces  an  acknowledged  error,  however,  we  may  be  (practically) 
unable  to  solve  the  original  non-linear  equations. 

3.  parameter  (or  statistical)  uncertainty:  provided  the  n  parameters  in  the  model  exist  and  are  unique, 
their  values  can  be  determined  from  n  experiments  if  both  the  model  and  experimental  data  are  exact. 
However,  in  the  presence  of  incomplete  information  (measurement  error,  model  error),  no  amount  of 
data  can  provide  exact  solutions  for  them  and  the  true  values  of  these  parameters  will  remain  uncertain. 

4.  Quantitative  Uncertainty  Modeling.  Uncertainty  quantification  can  be  accomplished  using  either  probabilis¬ 
tic  or  deterministic  methods  [28].  Examples  of  non-probabilistic,  sometimes  also  referred  to  as  possibilistic  methods, 
are  interval  analysis  [15],  Dempster-Shafer  theory  [36],  convex  modeling  [5],  and  fuzzy  computation  methods  [31]. 
Broadly  speaking  these  methods  address  the  possibility  of  an  event  within  the  bounds  imposed  on  the  uncertain  pa¬ 
rameters.  They  do,  however,  not  give  any  information  regarding  the  likelihood  of  this  event  taking  place. 

Probabilistic  methods,  on  the  other  hand,  provide  ample  quantitative  information  regarding  the  likelihood  of  an 
event  taking  place.  For  varying  degrees  of  computational  effort,  one  can  compute  statistics,  and  the  conhdence  inter¬ 
vals  thereof,  or  even  complete  probability  density  functions  (PDF)  for  any  particular  quantity  of  interest.  It  is  however, 
important  to  realize  that  the  model  outcomes  can  be  very  sensitive  to  the  accuracy  of  the  provided  probabilistic  in¬ 
put  [11].  It  is  important  to  recognize  that  not  just  marginal  PDFs  are  required  but  that  an  accurate  estimation  of  the 
correlation  between  random  variables  is  paramount  [41]. 

Most  specialists  agree  that  perhaps  the  most  promising  applications  for  the  non-probabilistic  methods  will  occur 
when  reliable  data  are  unavailable  or  hard  to  get.  Non-probabilistic  methods  are  mostly  used  in  absence  of  reliable 
and  accurate  data.  There  is  little  or  no  debate  that,  when  sufficient  data  exist,  the  probabilistic  model  is  superior. 
Probabilistic  computations  give  much  more  comprehensive  knowledge  about  the  state  of  a  particular  system.  The 
debate  between  the  believers  in  probabilistic  and  possibilistic  methods  will  probably  (or  is  it:  possibly?)  continue  for 
considerable  time.  In  short,  the  end-user  should  never  forget  that  no  matter  what  model  (probabilistic  or  possibilistic) 
is  used,  the  computed  outcomes  are  conditional  upon  the  model  assumptions  made  [12].  Uncertainty  associated  with 
the  models  should  not  be  taken  as  an  excuse  for  sloppiness!  No  matter  what  quantitative  uncertainty  analysis  model  is 
used,  the  validity  of  the  results  depends  on  the  validity  of  the  model. 

Recently,  a  mathematical  framework  which  unifies  the  probabilistic  and  possibilistic  algorithms  has  been  sug¬ 
gested  [19];  it  consists  of  a  general  constrained  optimization  problem: 

(4.1)  min{s(x)|(7(x)  =  0} 

where  g{x)  represents  the  limit  state  function,  i.e.  the  zero  safety  margin,  as  a  function  of  the  uncertain  parameters  x 
and  s(x)  is  some  function  which  characterizes  the  method.  For  instance,  when  s(x)  is  set  equal  to  the  reliability  index 
/),  the  solution  to  Eq.  (4.1)  reverts  to  the  well-known  first-order  reliability  method  (FORM)  [24]. 

In  this  paper,  only  probabilistic  analysis  methods  are  used.  Based  on  the  previous  discussion  of  uncertainty  types, 
the  following  3  issues  need  to  be  addressed  in  any  quantitative  uncertainty  assessment: 

1 .  characterize  uncertainty  associated  with  the  system  parameters  and  outside  environment 

2.  propagate  these  uncertainties  through  the  (potentially  large)  system 
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3.  incorporate  the  model  uncertainty  -lack  of  data  or  lack  of  knowledge  about  the  system  behavior-  into  the 
overall  assessment 

4.1.  Characterizing  aleatory  uncertainty.  A  complete  statistical  description  requires  the  identification  of  all 
marginal  and  joint  probability  density  functions  of  each  random  variable  in  the  problem  [13].  In  practice,  the  required 
amount  of  data  simply  will  not  be  available  or  the  data  structure  will  be  so  complex  that  only  an  approximate  descrip¬ 
tion  is  mathematically  tractable.  One  form  of  approximate  modeling  is  the  use  of  conditional  distributions,  where  the 
less  important  random  variables  are  modeled  conditional  upon  the  values  of  the  more  important  random  variables. 
Whether  a  particular  random  variable  is  important  depends  on  both  its  own  variability  and  on  the  sensitivity  derivative. 
Alternatively  a  second-order  description  can  be  used,  which  consists  of  the  marginal  PDF  for  each  of  the  variables 
and  the  correlation  coefficients  between  them.  The  correlations  can  be  hard  to  estimate  accurately  but  the  importance 
thereof  can  not  be  understated  and  is  clearly  illustrated,  for  instance,  by  the  expression  for  the  first-order  estimate  of 
the  variance  of  z  =  f{x,y) 


(4.2)  Var{*)  =  Var{x)  (g)  (g)  (|^) 

where  all  derivatives  are  evaluated  at  the  mean  values  lE(a;)  and  lE(y)  and  where 

(4.3)  Covar(a;,y)  =  VVar(a;)Var(y)  =  E  -  lE(a;)lE(y)] 


is  the  covariance  of  x  and  y  and  p^y  is  the  correlation  coefficient  of  x  and  y. 

In  the  absence  of  a  sufficient  amount  of  data,  expert  opinion  can  be  used  to  construct  an  approximate  joint  PDF.  An 
overview  of  expert  elicitation  techniques  is  described  in  [26]  and  [14].  It  should  be  stressed  however,  that  significant 
statistical  uncertainty  is  associated  with  such  techniques  (discussed  in  Section  4.3). 

4.2.  Propagating  uncertainties.  Consider  the  vector  x  =  {xi,. . .  ,x„)  and  its  joint  PDF  fx(xi,. ..  ,x„).  The 
relationship  between  the  input  variables  x  and  the  output  variables  y  =  {yi, ,  yn)  is  defined  by 


(4.4) 

together  with  the  inverse  functions 

(4.5)  a 
The  Jacobian  J  of  the  transformation  is: 

(4.6) 


Vi  =  gi{xi,...,x„)  for  i  = 


..-1 


J  = 


■  5  yn) 

for  i  =  f, 

dyi 

dyi 

dyi 

dxi 

dX2 

dXn 

dv2 

dy2 

dy2 

dxi 

dX2 

dXn 

dyn 

dyn 

dyn 

dxi 

dX2 

dXn 

The  exact  joint  PDF  for  the  model  output  y  =  g{x)  in  terms  of  the  y-variables  is  given  by  [16]: 


(4.7) 


fY{yi,---,yn) 


=  \J\  ^  fx{xi,...,Xn) 

=  \J\~^  fxioi^iyi,  -  ■  ■  ,yn),  ■  ■  ■  ,9n\yY  -  ■  ■  ,yn)) 


where  /,(•)  is  the  PDF  of  the  variable  •. 

Various  techniques  of  uncertainty  propagation  are  briefly  outlined  in  this  section.  Three  levels  of  methods  can  be 
identified: 


4 


1.  Simulation  methods:  In  their  most  elementary  form  (brute  force  Monte  Carlo)  one  samples  repeatedly  from 
the  joint  PDF  for  the  uncertain  variables  and  generates  the  histograms  of  the  model  response  quantities  of 
interest  [25] .  The  method  has  the  advantage  that  it  is  simple  and  universally  applicable.  The  convergence  does 
not  directly  depend  on  the  number  of  random  variables  in  the  problem.  In  this  form  Monte  Carlo  methods 
always  give  the  correct  answer,  but  a  prohibitively  large  number  of  simulations  may  be  required  to  accurately 
estimate  extreme  responses,  which  have  a  small  probability  of  occurrence.  Quite  often  these  are  the  ones 
of  interest  in  structural  engineering  and  for  this  reason  alternative  sampling  techniques  which  lead  to  much 
faster  convergence  have  been  developed.  An  overview  of  such  variance-reduction  techniques  is  available  in 
[13],  [25]  and  [34].  Such  techniques  should  be  used  carefully  since  importance  sampling,  e.g.,  can  introduce 
severe  bias  if  used  inappropriately. 

However,  the  brute  force  Monte  Carlo  method  requires  a  much  smaller  number  of  samples  -  typically  around 
100  or  less  -  to  accurately  estimate  the  mean  and  standard  deviation,  which  are  the  quantities  of  interest  in 
computational  fluid  dynamics  applications.  Brute  force  Monte  Carlo  also  sounds  appealing  since  its  imple¬ 
mentation  can  be  written  as  a  wrapper  around  an  existing,  deterministic  algorithm.  This  makes  it  easy  to 
transform  any  deterministic  algorithm  into  a  probabilistic  code.  However,  this  approach  is  not  very  efficient, 
particularly  when  the  deterministic  code  uses  an  iterative  solution  process.  The  combination  of  a  limited 
amount  of  extra  coding  to  make  use  of  restarts  and  a  smart  sampling  scheme  may  result  in  large  pay-offs  in 
execution  time. 

2.  Second-moment  analytic  approximations:  these  methods  are  based  on  a  Taylor  series  expansion  around  the 
mean  value  of  the  input  variables.  They  work  best  when  the  input  uncertainties  are  not  too  large  (coefficient 
of  variation  less  than  10  to  15%)  and  approximately  Gaussian.  The  methods  require  accurate  derivative 
information.  They  have  proven  highly  accurate  in  structural  mechanics  applications  [4]. 

3.  Advanced  analytic  methods:  These  methods  can  be  regarded  as  a  generalized  Monte  Carlo  method  [34]. 
However,  since  the  solution  mechanics  are  so  different,  we  treat  them  as  a  separate  category.  In  such  a 
method  the  stochastic  output  variables  are  written  as  a  series  expansion.  This  series  expansion  formulation 
readily  accommodates  advanced  random  held  modeling  of  the  uncertainties.  The  basic  equations  (ODE  or 
PDE)  are  now  solved  stochastically  and  a  large  amount  of  coding  is  required.  The  trade-off  is  that  an  accurate 
solution  can  be  obtained  for  moments  or  output  PDEs  which  accounts  for  both  the  non-linearities  in  the 
differential  equations  and  the  non-Gaussian  nature  of  some  of  the  input  variables.  Polynomial  Chaos  [39]  is 
one  member  of  this  class  of  methods. 

4.3.  Model  uncertainty.  Model  uncertainty  can  be  attributed  to  either  lack  of  data  or  their  inaccuracy,  or  an 
inexact  mathematical  model  form.  This  inexact  form  can  be  due  to  either  simplification  or  ignorance  [6].  Model 
parameters  cannot  be  estimated  accurately  when  insufficient  high-quality  data  are  available.  This  uncertainty  on 
the  parameters  can  be  reduced  by  collecting  more  data.  Due  to  the  asymptotic  normality  of  the  maximum  likelihood 
estimator,  a  second-order  description  or  confidence  bounds  of  the  statistical  uncertainty  is  readily  obtained  by  inverting 
the  local  Eisher  information  matrix  [18]. 

Model  inexactness  can  be  introduced  when  a  complex  model  is  simplified  and  some  of  the  variables  are  omitted 
from  the  analysis.  Sometimes  a  nonlinear  model  is  linearized  or  replaced  by  a  response-surface  model  to  save  computer 
time.  Uncertainties  due  to  ignorance  are  by  definition  notoriously  hard  to  quantify.  They  can  only  be  described  in 
probabilistic  terms  [12].  In  this  paper  we  will  consider  uncertainty  due  to  missing  variables  and  boundary  conditions. 

In  comparison  with  Eq.(4.7),  model  uncertainty  adds  an  extra  layer  of  uncertainty  to  the  problem,  but  does  not 
fundamentally  alter  the  structure  of  the  problem.  Let  6  and  M  denote  the  model  parameter  and  model  form  uncertainty, 
respectively.  When  the  model  uncertainties  6  and  M  are  not  taken  into  consideration,  the  results  obtained  using 
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inherent  variability  x  only  can  be  considered  conditional  upon  both  the  values  used  for  the  model  parameters  and  the 
mathematical  form  used  for  the  actual  model  itself:  y\6,M  =  g{x\6,M).  The  complete  PDF  model  for  the  output 
y  =  g{x,0,M)  is: 

= M) /e|M(^|M) /m(M) 

If  the  model  uncertainties  are  independent  of  the  inherent  uncertainties  the  complete  PDF  model  for  the  output  y  = 
g{x,  6,  M)  becomes  separable  and  simplifies  to: 

(4.9)  fY,s,Miy^S,M)  =  fY(y)fe(0)fM(M) 


5.  Burgers’s  Equation.  In  this  paper  we  consider  a  non-linear  1-D  model  problem,  which  includes  both  ad- 
vection  and  diffusion.  The  viscosity  is  the  only  model  parameter  in  this  problem.  The  complete  nonlinear  Burgers 
equation 


(5.1) 


du  du  _  d  /  du'\ 


is  a  parabolic  PDF  which  can  serve  as  a  model  equation  for  the  boundary-layer  equations  and  is  very  similar  to  the 
equations  governing  fluid  flow  [3].  For  simplicity  the  linearized  equation 


(5.2) 


du  du  _  d  f 


du  ,  ,  sdu  d‘^u 


is  often  used.  Equations  (5.1)  and  (5.2)  can  be  combined  into  a  generalized  equation  (GBE) 

(5.3) 

Eor  constant  viscosity  fi,  c  =  ^  and  b  =  —\,  the  GBE  (5.3)  has  the  stationary  solution 

(5.4)  u(a;)  =  i  [1  -b  tanh  j  j 


We  solved  the  GBE  as  a  steady-state  boundary  value  problem  with  a  second-order  accurate  central  difference  scheme 
using  a  conservative  formulation  for  Eq.(5.3).  The  non-linear  equation: 


(5.5) 


df  _  d  / 

d^^d^  V  ^ ) 


where  the  flux  /  =  (1  —  u),  is  solved  iteratively.  We  drive  the  residual  r{u)  to  zero: 


(5.6) 


df  d  (  du\ 


(5.7) 


dx  dx 

We  implemented  Newton’s  method  on  a  uniform  grid: 

(1^)^”^  Am(”)  = 

^(n+l)  ^  ^(n)  4.  ^^(n) 


Using  a  second-order  central  difference  scheme  this  leads  to  a  tri-diagonal  system  of  equations  in  each  iteration: 


hi  Cl 

{  K  > 

Aul 

r(ui) 

02  b2  C2 

0 

Au2 

r{u2) 

(5.8) 

0 

O-n  — 1  ^n  — 1  ^n  —  l 

< 

1 

<1 

>  =  -  < 

: 

Au„  ^ 

.  r{un)  , 
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Second  moment  estimates 


where 


(5.9) 


„  ,  _  9rj  _  (1/2-Uj-i) _ ^ 

duj-1  2Ax  Ax 

u  .  _  9^  _  M3+M3  +  1 
duj  A  3;^ 

_  dvj  _  (1/2-333  +  1)  _  M3  +  1 

3  9i3j  +  i  2A3; 


and  61  =  6„  =  1  and  ai  =  ci  =  a„  =  c„  =  r(ui)  =  r(un)  =  0. 


exact  deterministic  (ED) 


Figure  a 


First  Derivative 


Exact 


Fig.  6.1.  Comparison  of  exact  and  approximate  analytical  solutions  to  first  and  second  moment  ofu(x)  using  129  grid  points. 


6.  Random  Variable  Modeling  of  the  GBE.  The  differential  equation  Eq.(5.3)  is  solved  for  an  assumed  value 
of  the  viscosity  fi,  which  does  not  depend  on  the  independent  variable  x,  i.e.  ft  is  constant  over  the  entire  domain.  The 
exact  PDF  of  u{x)  can  be  obtained  from  the  following  equation  [16]: 


(6.1) 


fu{u{x))  = 


■^u{x,ii) 


-1 


fM{p) 


from  which  the  exact  mean  lE[u(a;)]  and  variance  Var[u(a;)]  =  lE[u^(a;)]  —  lE[u(a;)]^  are  readily  computed  using  Math- 
ematica  [42].  These  results  are  shown  in  Figure  6.1  for  an  assumed  a  Gaussian  distribution  for  ft  with  mean  value 
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E(/i)  =-p=  0.25  and  COV  of  10%. 

On  the  basis  of  the  closed-form  analytic  solution  in  Eq.(5.4),  we  can  compute  analytic  expressions  for  the  first 
and  second  derivatives.  We  can  readily  compare  first  and  second-order  approximations  with  exact  results.  The  first 
and  second  order  approximation  to  the  expected  value  E[u(a;)]  are  given  by: 


lEiro[M(a;)]  =  u{ji,x) 

E5:o[w(a;)]  =  u{jl,x)  +  5Var(/r)  _ 

The  first-order  approximation  is  identical  to  the  deterministic  solution.  The  mean  stochastic  and  second-order 
approximation  thereof  are  shown  on  Figure  6.1a.  The  graphs  almost  coincide  but  the  differences  E[u(a;)]  —  ^po  [w(a;)] 
and  E[u(a;)]  —  E5;o[w(a;)]  are  plotted  along  the  Y-axis  on  the  right  side  of  Figure  6.1a.  Since  the  error  between  the 
deterministic  (a.k.a.  first-order)  and  mean  stochastic  solution  has  a  shape,  which  is  very  similar  to  the  second-derivative 
(plotted  in  Figure  6.1b),  it  is  to  be  expected  that  the  second-order  correction  will  lead  to  a  much  more  accurate  result. 
The  error  on  the  second-order  first-moment  (SOFM)  estimate  is  2  orders  of  magnitude  smaller  than  the  error  on  the 
first-order  first-moment  (FOFM)  estimate.  It  is  worth  noting  that  it  follows  from  Fq.  (6.2)  that,  due  to  the  symmetric 
nature  of  the  Gaussian  PDF,  the  derivatives  of  odd-order  have  no  impact  on  the  estimate  for  E[u(a;)]. 

The  first  order  approximation  to  the  variance  of  u{x)  is: 


(6.3) 


ViiVFo[u{x)] 


When  /r  is  a  Gaussian  random  variable,  the  second  order  approximation  to  the  variance  of  u  is: 


(6.4) 


Var  so[u{x)] 


Var(/r)  -b  i 
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Both  the  exact  variance  Var[u(a;)]  and  the  FO  (6.3)  and  SO  (6.4)  approximations  are  even  functions  of  x.  Figure 
6. Id  compares  the  accuracy  of  the  first-order  and  second-order  solutions.  It  can  be  concluded  that  the  second-order 
estimate,  which  is  always  greater  than  the  first-order  estimate  for  a  Gaussian  random  variable,  is  generally  speaking 
more  accurate.  Figure  6.2  shows  the  true  PDF  of  u{x)  for  selected  values  of  x.  The  PDF  shapes  for  u{xo)  and  u{—xo) 
are  identical,  but  mirrored.  Both  the  FOSM  and  SOSM  estimate  assume  that  the  underlying  distribution  is  symmetric. 
Consequently,  the  largest  errors  occur  near  the  ends  of  the  a;-interval,  i.e.  where  the  PDF  for  u{x)  is  strongly  skewed. 
Since  the  PDF  at  a;  =  ±1  is  nearly  Gaussian,  the  SOSM  estimate  is  highly  accurate.  It  is  also  important  to  note  that 
a  correct  value  Var[u(0)]  =  0  is  predicted  by  all  methods.  The  skewness  and  kurtosis  are  shown  in  Figure  6.3.  The 
kurtosis  value  /t  =  0  at  a;  =  0  is  caused  by  the  degenerate  PDF  for  u(0). 


7.  Gaussian  Random  Fields. 


7.1.  Second-Order  Description.  A  continuous  (or  discrete)  random  field  can  be  regarded  as  an  infinite  (or 
finite)  set  of  random  variables  defined  over  its  domain  X.  One-dimensional  random  fields  Z{x)  are  also  referred  to 
as  random  processes.  For  all  practical  purposes,  continuous  random  fields  must  be  discretized  for  use  in  numerical 
algorithms,  and  we  can  concentrate  on  discrete  random  field  characteristics  only.  Just  like  any  other  n-dimensional 
probability  density,  a  random  field  is  fully  characterized  either  by  its  joint  density  or  by  its  marginals  and  all  its  higher- 
order  conditional  PDFs,  i.e  of  order  2, 3, . . . ,  n  —  1  [16].  In  second-order  analysis,  a  random  field  is  described  by  the 
marginals  and  the  auto-correlation  function  R{x,  x'),  which  is  usually  defined  as  the  correlation  between  the  random 
variables  Z{x)  and  Z{x')  at  two  locations  x  and  x'  [40]: 


,  ^  Covar[Z{x),Z{x')] 
^Vax[Z{x)]-Vax[Z{x')] 


(7.1) 


probability  density  function  probability  density  function 


Figure  c:  x  =  +3  Figure  d:  a;  =  +1 

Fig.  6.2.  Exact  PDF  of  u(x)  at  selected  values  for  x. 


For  a  homogeneous  random  field,  the  marginal  distribution  is  independent  of  the  location  and  the  auto-correlation 
function  R  depends  on  only  the  distance  or  lag  ^  =  x  —  x'  between  the  two  locations  [29].  A  homogeneous,  Gaussian 
process  is  completely  determined  by  its  mean  value  lE[Z(a;)],  variance  Var[Z(a;)]  and  auto-correlation  R{^). 

Table  7.1 

Commonly  used  families  of  auto-correlation  models  and  corresponding  spectral  density  junctions  [29]. 


Correlation  model 

Auto-correlation  R{^) 

Spectral  density  S{u;) 

exponential 

exp(-|^|//c) 

lo 

squared  exponential 

exp  [- 

2^exp[-(^)"] 

triangular 

i  1-S  if 

1— cos(c<;/c) 

~kuPT~c 

o 

1— h 

V 

(  k-  if  |w|  <  f 

damped  sinusoidal 

sin(  7r\^\  /Ici 

Tr\^\/lc 

1  27t  I^I  — 

1  0  if  |c.|  >  ^ 

Figure  a:  compare  with  71  =  0  for  Gaussian  PDF  Figure  b:  compare  with  /t  =  3  for  Gaussian  PDF 

Fig.  6.3.  Skewness  and  kurtosis  of  the  exact  PDF  ofu(x)  using  129  grid  points. 


Some  commonly  used  models  for  the  auto-correlation  function  are  listed  in  Table  7.1  and  shown  in  Figure  8.2a. 
The  first  two  are  perhaps  most  commonly  encountered  in  the  literature  [21].  The  first  one  applies  to  a  non-differentiable 
field,  whereas  the  second  one  is  found  in  differentiable  fields.  The  correlation  length  1^  is  a  very  important  measure  for 
the  rate  of  fluctuation  within  the  random  held  model.  A  larger  correlation  length  1^  corresponds  to  a  slower  varying 
random  held  as  is  illustrated  in  Figure  7.1  for  the  exponential  correlation  model. 


Fig.  7.1.  Exponential  auto-correlation  function  and  sample  random  fields  for  different  correlation  length  Ic-  Marginal  density  for  the  viscosity 
is  Gaussian  with  mean  value  0.25  and  COV=  10%. 


7.2.  Discretization  of  Random  Fields.  Most  modem  numerical  solution  algorithms  to  PDF’s  (such  as  finite 
element,  finite  volume  or  finite  difference)  require  that  all  continuous  parameter  fields  are  discretized.  When  spatial 
or  temporal  uncertainty  effects  are  directly  included  in  the  analysis  it  makes  sense  to  use  random  fields  for  a  more 
accurate  representation  of  the  variations.  For  instance,  a  thickness  of  the  wing-skin  is  random  but  is  also  not  likely  to 
be  exactly  the  same  everywhere.  The  spatial  or  temporal  fluctuations  of  a  parameter  cannot  be  accounted  for  if  that 
parameter  is  modeled  by  only  a  single  random  variable. 

We  only  give  a  brief  overview  of  some  of  the  discretization  methods,  a  more  extensive  review  is  available  in  the 
literature  [34].  Some  of  these  methods  are  limited  to  Gaussian  fields  only.  We  will  discuss: 
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1 .  midpoint  method 

2.  locally  averaged  method 

3.  shape  function  method 

4.  series  expansion  method 

Other  methods  include  locally  averaged  subdivision  [10]  and  the  turning-bands  method  [9]. 

The  simplest  method  of  discretization  is  the  midpoint  method  [7].  In  this  method  the  held  within  the  one  element 
or  cell  is  described  by  a  single  random  variable,  which  represents  the  held  value  at  the  centroid  of  the  element.  The 
value  of  the  held  is  assumed  to  be  constant  within  the  entire  cell.  Consequently,  the  discretized  held  is  step-wise 
constant  with  discontinuities  at  the  cell  boundaries.  Its  statistics  (mean,  variance,  . . .)  are  readily  computed  from  the 
statistics  at  the  centroids. 

In  the  locally  averaged  method,  the  held  within  each  cell  is  described  in  terms  of  the  spatial  average  of  the  held 
over  the  element  [40].  The  discretized  held  is  still  constant  over  each  cell  with  step-wise  discontinuities  at  the  cell 
boundaries,  but  one  can  expect  a  better  ht  due  to  the  averaging  process.  Typically,  the  variance  of  the  averaged  random 
process  is  smaller  than  the  local  variance  of  the  random  held,  used  for  midpoint  discretization.  Expressions  for  the 
statistics  will  be  derived  in  Section  9. 

The  shape  function  method  [23]  describes  the  random  held  using  shape  functions  anchored  at  a  set  of  nodal  values. 
The  method  is  particularly  elegant  and  efficient  in  conjunction  with  the  hnite  element  method  [22].  The  random  held 
description  is  continuous  across  cell  boundaries  and  this  represents  a  clear  improvement  over  the  midpoint  and  locally 
averaged  discretization  methods. 

Yet  another  method  uses  a  series  expansion  approach  to  discretize  the  random  held.  Legendre  basis  functions  are 
used  with  Gaussian  random  helds  in  [20].  Other  functions  can  be  used  in  conjunction  with  different  marginal  PDFs 
[43].  Provided  that  the  exact  eigenfunctions  in  the  Karhunen-Loeve  expansion  are  available,  this  method  requires  the 
smallest  number  of  random  variables  for  a  given  level  of  accuracy  [39].  However,  often  there  is  no  exact  solution  to  the 
integral  eigenvalue  problem  associated  with  the  Karhunen-Loeve  expansion.  In  this  case,  the  resulting  approximations 
essentially  reduce  the  series  expansion  method  to  a  shape  function  method  [21]. 

In  this  report,  we  compare  results  using  the  midpoint  and  locally  averaged  random  held  discretization  methods 
for  Monte  Carlo  simulations.  Polynomial  chaos  results  will  be  reported  in  an  upcoming  report  [41]. 


8.  Generation  of  Midpoint  Random  Field  Samples.  Generally  speaking,  three  families  of  methods  can  be  used 
to  generate  sample  random  helds.  The  different  methods  will  be  applied  only  in  conjunction  with  the  auto-correlation 
functions  given  in  Table  7.1,  but  are  -  in  most  cases  -  applicable  to  other  auto-correlation  functions  as  well.  The 
random  held  can  be  discretized  over  the  grid,  which  is  used  to  solve  the  differential  equation,  or  over  any  other  grid. 
Typically,  the  random  held  grid  is  a  subset  of  the  hnite  element/difference/volume  grid  but  this  -  strictly  speaking  - 
need  not  be  the  case.  Only  Gaussian  random  helds  are  included  in  this  discussion. 


8.1.  Auto-regressive  processes.  Processes  with  an  exponentially  decaying  auto-correlation  function  can  easily 
be  transformed  into  hrst-order  auto-regressive,  or  AR(1)  processes,  which  are  also  known  as  Markov  processes.  In 
such  processes  the  auto-correlation  decreases  steadily  with  increasing  distance  between  sample  points.  For  a  uniform 
grid  the  discrete  AR(1)  process  Z{x)  is  dehned  by: 

(8.1)  Z{xi)  =4)Z{xi-i)+ai 
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correlation  coefficient 


where  |^|  <  1  and  all  a,  are  independent,  identically  distributed  zero-mean  Gaussian  variables.  The  discrete  process 
has  the  following  second-order  characteristics  [35]: 


^[Z{x)]  =  0 

(8.2)  y^[Z{x)]  =  ^ 

Rz{xi,Xi±k)  = 

If  we  chose  the  parameters 

(8  3)  §!-  =  exp(-l/g 

Var[aj]  =  1  -  exp(-2/(c) 

it  readily  follows  that  the  field  Z{x)  has  an  exponential  correlation  function  Rz{xi,Xi±k)  =  Rz{0  —  6xp(— |^|/(c), 
where  ^  =  Xi  -Xi^u- 


Figure  a:  correlation  function  Figure  b:  L2  Norm  of  error  estimate 

Fig.  8.1.  Theoretical  and  sample  auto-correlation  functions  for  exponential  correlation  modelwithlc  =  1. 


Table  8.1 

Mean  Relative  (in  percent)  and  Absolute  Squared  Error  on  the  sample  statistics  for  the  exponential  correlation  model  with  Ic  =  1.  Marginal 
density  for  the  viscosity  is  Gaussian  with  mean  value  0.25  and  COV=  10%. 


^samples 

Relative  error  (in  %)  on 

Mean  St.  Dev.  Correl. 

Absolute  error  (L2  norm)  on 

Mean  St.  Dev.  Correl. 

100 

0.14 

0.61 

1.13 

0.00035 

0.00015 

0.00503 

1000 

0.02 

0.24 

1.43 

0.00006 

0.00006 

0.00605 

10000 

0.01 

0.07 

0.06 

0.00002 

0.00002 

0.00029 

For  uniform  grids,  the  sample  auto-regressive  processes  are  easily  generated  using  this  method.  The  convergence 
of  the  sample  statistics  to  the  true  values  of  the  random  field  is  shown  in  Figure  8.1  and  Table  8.1.  From  a  practical 
point  of  view,  the  stability  and  accuracy  of  the  ensemble  statistics  is  very  important  for  the  simulation  process.  A  more 
general  description  of  auto-regressive  processes  can  be  found  in  [33]  and  [27]. 
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8.2.  Covariance  decomposition.  In  this  method,  the  covariance  matrix  of  the  discretized  field  is  generated  first. 
For  a  zero-mean,  homogeneous  process  Z  (x)  the  symmetric,  positive  definite  covariance  matrix  is  directly  derived 
from  the  auto-correlation  function: 

(8.4)  Cov£iv[Z{xi),Z{xj)]  =  i?(^)Var(Z) 

where  ^  =  ja;,  —  |  is  the  distance  between  the  points  Xi  and  xj. 

Along  with  the  mean  value  lE[Z(a;)]  =  Z,  the  covariance  matrix  fully  describes  the  discretized  Gaussian  random 
field.  Using  the  covariance  decomposition  method  the  sample  values  Zj{xi)  of  Z{xi),  i  =  1 ..  .n  can  be  generated  on 
the  basis  of  n  independent  random  normal  variables  yc 

(8.5)  {Z{xi)}  =  {Z]  +  [A\^.{yi} 

where  the  matrix  [A]  contains  the  eigenvectors  of  the  covariance  matrix.  The  variance  of  each  of  the  independent 
Gaussian  variables  y,  is  equal  to  A,,  the  eigenvalue  of  the  covariance  matrix.  In  this  method  the  sample  mean  and 
covariances  will  only  approximately  match  the  required  values  and  approach  the  true  values  as  the  number  of  samples 
increases. 

It  is  possible  to  ensure  that  each  ensemble  matches  the  prescribed  means  and  covariances,  even  with  a  small 
number  of  samples.  This  preconditioning  method  [44]  combines  the  modal  decomposition  of  the  covariance  matrix 
with  the  spectral  representation  of  the  random  field,  which  is  described  in  the  next  section.  However,  in  this  method 
the  variates  y,  or  Z{xi)  are  not  directly  generated  from  a  Gaussian  distribution;  they  are  only  approximately  Gaussian 
by  means  of  the  central  limit  theorem.  They  are  generated  as  a  sum  of  independent  and  identically  distributed  random 
variables,  which  approaches  a  Gaussian  random  variable  if  sufficient  terms  are  included  in  the  sum. 

The  covariance  decomposition  method  readily  identifies  the  most  important  random  variables  as  the  ones  with 
the  largest  eigenvalues.  This  is  useful  if  one  wants  to  reduce  the  number  of  random  variables  in  the  random  field 
description. 

8.3.  Spectral  representation.  Unlike  the  previous  methods,  the  representation  of  a  random  field  using  the  spec¬ 
tral  method  is  continuous  across  the  elements.  The  underlying  idea  is  that  a  continuous  random  function  can  be 
expanded  on  a  basis  of  deterministic  functions  with  random  coefficients.  The  goal  is  to  achieve  a  sufficiently  accurate 
representation  of  the  random  field  with  a  finite,  i.e.  discrete,  number  of  functions  and  coefficients. 

The  spectral  density  S(lo)  (or  power  spectrum)  of  a  stationary  process  Z{x)  is  given  by  the  Fourier  transform  of 
its  auto-correlation 

(8.6)  5(w)  =  —  /  i?(^)exp(-jw^)d^ 

J  —  OO 

Since  i?(^)  is  an  even  function  for  a  real-valued  process  Z{x),  it  follows  that  the  S{uj)  is  a  real  function  of  cu: 

1 

(8.7)  S(lo)  =  -  /  R(0  cosiioOd^ 

^  Jo 

The  random  field  discretization  is  achieved  through  the  discretization  of  the  spectral  density.  The  power  spectrum 
is  approximated  at  some  pre-selected  discrete  frequencies,  typically  chosen  equally  spaced  in  [— Wmax,  i^max]- 

_  Nf,,, 

(8.8)  Zj{x)  =  Z  +  ^2  'J‘Z^^^[Z{x)]S{ijJi)/S.ijjcos{ijj\x  +  (fy) 

where  Aw  is  the  spacing  between  the  sampling  frequencies,  w'  is  a  small  perturbation  of  w,  and  is  the  phase  angle 
[37]. 
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The  process  generated  according  to  Eq.(8.8)  is  asymptotically  Gaussian  as  Nf^eq  oo.  Since  the  highest  fre¬ 
quencies  are  ignored  in  the  discretization  of  the  random  field,  part  of  the  variance  of  the  process  will  be  missed. 
Another  source  of  discretization  error  stems  from  the  finite  number  of  terms  used  in  the  representation.  The  conver¬ 
gence  to  the  target  auto-correlation  R{0  or  to  the  target  spectral  density  S{u;)  is  inversely  proportional  to  [38]. 

Additional  details  of  this  method  are  described  in  [37]. 

If  the  eigenvalues  and  eigenfunctions  of  the  covariance  matrix  Eq.(8.4)  are  known,  the  Karhunen-Loeve  expan¬ 
sion  requires  the  smallest  number  of  functions  to  achieve  a  prescribed  level  of  accuracy.  However,  in  many  cases 
the  eigenfunctions  are  only  approximately  known  and  the  approximate  K-L  expansion  is  no  more  efficient  than  the 
shape  function  expansion  [22].  Other  expansions  may  require  more  terms  to  achieve  the  same  accuracy  but  require 
substantially  less  computational  effort.  A  detailed  discussion  of  the  truncation  error  in  each  of  the  methods  can  be 
found  in  [21]  and  [45]. 

The  slow  fluctuations  are  represented  using  the  lowest  frequencies,  whereas  the  high  frequencies  are  required  to 
accurately  model  the  rapid  fluctuations.  It  automatically  follows  that  much  higher  frequencies  are  required  in  Eq.(8.8) 
to  accurately  model  the  sample  random  held  in  Eigure  7.1a  than  for  the  one  in  Eigure  7.1c.  Eigure  8.2  illustrates  the 
relationship  between  the  auto-correlation  model  and  the  frequency  range  which  is  required. 

9.  Locally  Averaged  Random  Fields.  Eor  the  1  -D  problem  considered  here  we  get  the  following  statistics  for 
the  locally  averaged  random  held; 

1  l‘X-\-Xj2 

(9.1)  Zx{x)  =  -  /  Z{s)ds 

^  Jx-X/2 

where  X  is  the  cell  length  and  x  is  the  midpoint  of  the  cell.  Since  the  held  is  homogeneous, 

(9.2)  W.{Zx{x))=W.{Z{x))=^ 

The  averaging  process  typically  reduces  the  variance  [40].  The  variance  function  7(A)  measures  the  reduction  of  the 
point  variance  Var[Z]  due  to  the  averaging  process.  Eor  a  homogeneous  random  held;; 

(9.3)  Var[Zx]  =  =  7(A)Var[Z] 


Eigure  a:  Auto-correlation  models 


Eigure  b:  Power  spectral  density 


Fig.  8.2.  Auto-correlation  functions  and  corresponding  spectral  densities  used  in  simulations. 
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and,  for  a  homogeneous  field,  the  variance  function  solely  depends  on  the  cell  size  X : 


(9.4) 


7(X) 


(^1  -  RiOd^ 


Fig.  9.1.  Definition  of  intervals  and  distances  used  in  Eq.(9.7). 

The  scale  of  fluctuation  6  is  defined  as  follows: 

(9.5)  e=  lim  Xj{X) 

X  — >-00 

which  indicates  that,  asymptotically  for  large  averaging  length  X:  ^{X)  ~  9  jX.  Combination  of  Eqs.(9.5)  and  (9.4) 
implies  that 

poo  poo 

(9.6)  9  =  2  R(Od^  =  /  R(\^\)d^ 

Jo  J  —  oo 

It  follows  immediately  from  the  spectral  density  definition  in  Eq.  (8.6)  that  the  scale  of  fluctuation  is  proportional  to 
the  ordinate  5(a;  =  0)  of  the  spectral  density  function  9  =  27r5(0). 

The  scale  of  fluctuation  provides  a  lot  of  information  over  the  extent  of  the  correlation  within  the  random  field.  If 
the  correlation  dies  out  rapidly,  the  asymptotic  relationship  will  be  achieved  for  a  relatively  small  averaging  length  X. 
If  two  locally  averaged  field  values  are  more  than  distance  9  apart,  Eq.(9.5)  indicates  that  they  are  almost  independent 
of  each  other. 

On  the  basis  of  Eq.  (9.4),  the  second-moment  properties  of  the  locally  averaged  random  field  can  be  computed. 
These  can  then  be  used  to  generate  sample  random  fields  using  the  covariance  decomposition  method.  The  correlation 
coefficient  pz^  ,z^i  between  the  locally  averaged  random  field  values  Zx  and  Zx'  is  given  by: 


(9.7) 


PZxZx' 


A{x  -X)  +  A{x  +  X)-  2A{x) 
2A{X) 


where  A(«)  =  (•)^7(«)  and  the  distances  x  and  X  =  X'  sxe  defined  in  Eigure  9.1. 


10.  Random  Field  Modeling  of  the  GBE.  In  reality,  the  viscosity  also  depends  -  among  other  factors  -  on  the 
temperature,  for  which  no  information  is  available  in  the  Burgers  model  (5.3).  The  model  uncertainty  caused  by  these 
missing  variables  is  best  described  using  a  random  field,  which  can  be  thought  of  as  an  infinite  set  of  random  variables. 
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each  one  of  them  describing  the  variability  of  the  viscosity  at  a  particular  location  x.  The  second-order  description  of 
a  homogeneous  random  field  [29]  is  given  by  the  marginal  PDF  /(/t)  and  the  auto-correlation  function  R{^)  which 
are  both  independent  of  the  location  x  for  a  homogeneous  random  field: 


(10.1) 


^  Covar[/r(a;),/t(a;-b^)] 

yar[^l{x)] 


where  ^  =  X2  —  xi  is  the  lag  or  distance  between  the  two  locations  considered  [40].  It  is  assumed  that  the  marginal 
density  /(/t)  is  Gaussian  with  mean  value  Jt  =  0.25  and  a  COV  of  10%. 

The  impact  of  the  uncertain  boundary  conditions  on  the  solution  is  discussed  in  the  next  section.  This  uncertainty 
is  ignored  in  this  section:  we  directly  applied  the  Dirichlet  boundary  conditions,  obtained  from  the  deterministic 
steady-state  solution  Eq.(5.4): 

(10.2)  u(±a;max)  =  ^  [l  +tanh 


-  Deterministic  -  Random  Field 

—  —  —  -  Random  Variable  -  Random  Field 


X 


Fig.  10.1.  Difference  between  deterministic,  mean  random  variable  and  mean  random  field  solution  (exponential  auto-correlation  with  1^  =  1). 

The  results  in  Section  6  using  a  random  variable  model  reveal  that  the  difference  between  the  deterministic  and 
the  mean  stochastic  solution  is  small  from  a  practical  point  of  view.  This  is  still  the  case  for  the  random  field  model  as 
well.  Figure  10.1  indicates,  however,  that  the  difference  between  the  deterministic  and  the  mean  random  field  solution 
is  somewhat  smaller  than  the  difference  between  the  mean  random  variable  and  the  mean  random  field  solution.  To  a 
large  extent,  this  can  be  attributed  to  the  formulation  of  the  boundary  conditions  using  Eq.(10.2);  which  are  identical 
for  the  deterministic  and  random  field  solution.  Note  that  the  mean  random  variable  solution  is  an  exact  result,  whereas 
the  mean  random  field  solution  is  the  results  of  10®  Monte  Carlo  simulations.  The  standard  error  on  the  Monte  Carlo 
results  is  less  than  0.00006. 

Eor  this  reason  we  will  not  study  the  impact  of  the  correlation  length  1^  and  auto-correlation  model  i?(^)  on  the 
mean  solution.  The  discussion  is  limited  to  the  effects  of  the  random  field  modeling  on  the  standard  deviation  only. 
It  immediately  follows  from  Eigure  7.1  that  the  random  field  and  random  variable  models  lead  to  identical  solutions 
(shown  in  Eigure  6.1c)  for  an  infinite  correlation  length  1^  =  00. 

10.1.  Comparison  of  midpoint  and  locally  averaged  solutions.  When  the  viscosity  fj,{x)  is  represented  by  a 
random  field,  the  numerical  solution  of  the  GBE  in  Eq.(5.3)  requires  a  discretization  at  two  different  levels: 

1 .  discretization  required  to  solve  the  differential  equation  and  compute  an  accurate  solution  for  the  velocity 
field  u{x). 
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Standard  deviation  of  u  at  x  =  >0.75 


Viscocity  modeled  as  a  random  field  with  exponential  auto-correlation 
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—  0 —  4  s  0.001 
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1 0000  samples,  error  on  standard  deviation  is  less  than  3%  with  95%  confidence 
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Viscocity  modeled  as  a  random  field  with  exponential  auto-correlation 


Figure  a:a;  =  —0.75 


Figure  b:  a;  =  0 


Fig.  10.2.  Cell  size  dependence  of  the  variance  using  the  midpoint  discretization  method. 


2.  discretization  required  to  model  the  fluctuations  within  the  random  held  ij{x)  and  ensure  that  the  resulting 
discretized  equations  are  an  accurate  representation  of  the  stochastic  differential  equation. 

A  more  detailed  error  analysis  of  the  second-order  central  difference  scheme  used  in  the  analysis  is  documented 
elsewhere  [41].  In  this  paper,  the  focus  is  on  the  discussion  of  the  required  discretization  level  of  the  random  held 
itself.  All  the  results  in  this  Section  apply  to  the  exponential  correlation  model,  but  the  same  conclusions  can  be  drawn 
for  other  correlation  models. 

In  this  section  we  compare  Monte  Carlo  results  obtained  using  either  the  midpoint  or  the  locally  averaged  dis¬ 
cretization  method.  The  sample  random  fields  for  the  midpoint  method  are  generated  using  the  auto-regressive  property 
of  the  exponential  correlation  model  Eq.(8.1).  The  sample  random  fields  for  the  locally  averaged  method  are  obtained 
using  covariance  decomposition  Eq.(8.5). 

It  is  intuitively  clear  [7]  that  the  midpoint  discretization  method  tends  to  overestimate  the  variability  of  the  output 
quantities  of  interest  if  the  discretization  is  too  coarse.  This  can  be  explained  as  follows.  When  the  cell  size  is  much 
larger  than  the  correlation  length  1^  or  the  scale  of  fluctuation  6,  the  viscosity  fluctuates  rapidly  with  the  cell  (like  in 
Eigure  7.1a).  Consequently,  the  midpoint  values  of  the  different  cells  will  be  practically  independent  of  each  other.  As 
a  result  thereof,  Var[u(a;)]  will  decrease  if  the  number  of  cells  increases.  Eor  uncorrelated  midpoint  values  Var[u(a;)] 
will  be  inversely  proportional  to  the  number  of  cells.  Eor  a  non-zero  correlation  length  1^,  the  fluctuations  virtually 
disappear  when  the  cell  size  is  smaller  than  the  correlation  length  1^  (like  in  Eigure  7.1c).  In  this  case  the  midpoint 
values  of  neighboring  cells  become  highly  correlated.  Dividing  the  cells  even  further  leads  to  basically  identical 
midpoint  values  for  ij{x)  in  the  newly  created  cells  and  Var[u(a;)]  stabilizes. 

Eigure  10.2  illustrates  this  effect  for  various  correlation  lengths  1^  and  an  exponential  correlation  model  for  two 
locations: 

•  X  =  —0.75:  this  is  near  the  location  where  the  random  variable  model  (Eigure  6.1c)  has  its  largest  variance. 

•  a;  =  0:  this  is  the  location  where  the  random  variable  model  (Eigure  6.1c)  has  zero  variance. 

When  11  is  modeled  as  a  random  variable,  the  Burgers  equation  is  solved  once  for  an  assumed  value  of  the  viscosity 
II.  The  mathematical  expression  of  the  closed  form  solution  for  u{x),  given  in  Eq.(5.4),  is  independent  of  the  actual 
value  of  II.  When  the  viscosity  is  modeled  as  a  random  held,  however,  the  ODE  in  Eq.(5.3)  becomes  a  stochastic 
differential  equation.  It  is  immediately  clear  from  Eigure  10.2b  that  the  solution  of  the  stochastic  differential  equation 
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Standard  deviation  of  u  at  x  =  >0.75 


has  entirely  different  properties  than  Eq.(5.4).  The  random  variable  solution  of  Eq.(5.4)  results  in  Var[u(0)]  =  0. 
Eigure  10.2b  indicates  that  this  result  is  recovered  in  the  limit  for  1^  =  oo,  but  that  this  is  generally  not  the  case  for  a 
finite  correlation  length. 

Eigure  10.2  shows  that  the  standard  deviations  stabilizes  once  the  cell  size  is  roughly  equal  to  the  scale  of  fluc¬ 
tuation  0.  Recall  that  9  =  21^  for  the  exponential  auto-correlation  model.  It  is  interesting  to  note  that  the  white  noise 
random  field  model  (1^  =  0)  converges  to  a  deterministic  response,  i.e.  Var[u(a;)]  =  0. 


Results  for  locally  averaged  and  midpoint  discretization  of  viscocity 
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10000  samples,  error  on  standard  deviation  is  less  than  3%  with  95%  confidence 
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Results  for  locally  averaged  and  midpoint  discretization  of  viscocity 
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10000  samples,  error  on  standard  deviation  is  less  than  3%  with  95%  confidence 
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Fig.  10,3.  Comparison  of  cell  size  dependence  of  the  variance  using  the  midpoint  and  the  locally  averaged  discretization  method. 

On  the  other  hand,  locally  averaged  discretization  causes  the  true  variance  to  be  underestimated  [40].  When  the 
cell  size  is  much  larger  than  1^  (as  in  Eigure  7.1a)  the  local  averaging  process  drastically  reduces  the  variance  as 
explained  in  Section  9  (variance  function  7  close  to  0).  Consequently,  the  variance  of  the  discretized  viscosity  field 
in  Eq.(5.3)  will  be  underestimated.  When  the  cell  size  is  much  smaller  than  1^  (as  in  Eigure  7.1a)  the  local  averaging 
process  leaves  the  variance  virtually  intact  (variance  function  7  close  to  1)  and  the  full  variance  of  the  process  is 
included  in  the  discretized  random  field.  Eigure  10.3  indicates  that  a  stable  solution  for  the  variance  is  obtained  if  the 
cell  size  is  smaller  than  the  correlation  length  1^  and  is  approached  from  below  as  the  cell  size  decreases.  Eigure  10.3 
also  shows  that  the  midpoint  and  locally  averaged  results  will  tend  to  bracket  the  true  variability  of  the  output  quantity 
of  interest. 

Eigure  10.4  summarizes  the  impact  of  the  correlation  length  1^  on  Var[u(a;)].  It  follows  from  the  figure  that  for 
practical  purposes  it  suffices  to  only  know  the  order  of  magnitude  of  l^.  Eor  1^  00,  the  random  field  is  represented 

by  a  single  random  variable  and  we  obtain  the  same  solution  as  in  Eigure  6.1c. 

10.2.  Effect  of  auto-correlation  function  or  spectral  density.  In  this  section,  the  effect  of  the  auto-correlation 
model  on  the  standard  deviation  is  assessed.  All  Monte  Carlo  solutions  are  computed  using  midpoint  discretizations 
of  the  viscosity  field.  The  grids  are  fine  enough  to  obtain  converged  values  for  the  statistics.  All  random  field  sam¬ 
ples  are  generated  using  the  spectral  representation  method.  The  power  spectral  densities  were  truncated  at  Wmax, 
corresponding  to  the  95%-level: 

pOJrnax 

(10.3)  2  /  S{Lo)dLo  =  0.95 

Jo 

Eor  all  computations  the  frequency  range  [— Wmaxji^max]  is  discretized  using  100  evenly  spaced  frequencies;  the 
truncation  frequency  for  each  model  is  listed  in  Table  10.1.  The  standard  deviations  computed  for  two  different 
assumed  correlation  lengths  =  0.1  and  1^  =  10  using  all  4  auto-correlation  models  are  shown  in  Eigure  10.5. 
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standard  deviation  of  u(x) 


1 0000  samples,  error  on  standard  deviation  is  less  than  3%  with  95%  confidence 


Fig. 


10.4.  Impact  of  the  correlation  length  Ic  on  the  standard  deviation 


(midpoint  discretization  method,  all  results  are  grid-converged) 


Figure  a:  =  0.1 


Fig.  10.5.  StDev[t((a;)]  obtained  using  1000  midpoint  Monte  Carlo  simulations  of  the  random  field  and  512  grid  cells. 


As  suggested  by  Figure  10.4  for  the  exponential  auto-correlation  model,  Figure  10.5  shows  that  for  all  correlation 
kernels  used  the  variability  is  generally  larger  for  1^  =  10  than  for  1^  =  0.1.  For  the  shorter  correlation  length  =  0.1 
the  triangular  and  damped  sinusoidal  correlation  model  lead  to  nearly  identical  variability  in  u{x)  (Figure  10.5a). 
Figure  10.5b  shows  that  this  is  by  no  means  a  universal  trend,  but  depends  on  the  correlation  length  l^.  It  can  be 

Table  10.1 

Truncation  frequency  Umax  used  for  computation  of  results  in  Figure  10.5. 


Correlation  model 

^max^c 

exponential 

12.7 

squared  exponential 

3.92 

triangular 

13.0 

damped  sinusoidal 

2.98 
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concluded  that  for  l^.  =  0.1,  the  fluctuations  in  the  viscosity  held  are  too  rapid  for  the  weaker  medium-to-long-range 
correlations,  caused  by  oscillating  behavior  of  the  correlation  function  (see  Figure  8.2a),  to  have  an  effect.  The 
weak  correlations  between  points  separated  by  a  distance  larger  than  1^  are  drowning  within  the  rapidly  fluctuating 
field.  When  the  viscosity  held  is  varying  more  slowly,  these  correlations  are  strong  enough  to  have  an  impact  on  the 
solution.  Note  that  the  slight  asymmetry  in  Figure  10.5b  (most  pronounced  for  the  squared  exponential  model)  is 
entirely  statistical  in  nature.  Only  1000  Monte  Carlo  simulations  were  used  to  compute  the  results;  the  asymmetry 
disappears  when  more  simulations  are  performed. 

11.  Uncertain  Boundary  Conditions.  In  the  previous  section,  deterministic  Dirichlet  boundary  conditions  were 
imposed.  In  many  CFD  applications  the  differential  equations  have  to  be  integrated  over  a  (semi)-inflnite  domain. 
In  practice,  the  integration  domain  is  cut  off  at  some  “far  held”  boundary  and  deterministic  boundary  conditions  are 
imposed.  This  approach  essentially  ignores  all  variability  outside  the  truncated  integration  domain.  In  this  section, 
the  impact  of  the  uncertainty  modeling  of  the  boundary  conditions  on  the  second-order  statistics  of  the  solution  in  the 
interior  of  the  integration  domain  is  analyzed. 

In  many  practical  cases,  an  exact,  statistical  treatment  of  the  effects  of  the  truncation  of  the  integration  domain 
will  be  impossible.  Also,  the  data  required  to  build  the  sophisticated  statistical  model  for  such  an  analysis  may  be 
unavailable.  Therefore,  we  focus  on  a  more  practical  second-order  treatment  of  these  uncertainties. 

In  particular  we  will  compare  the  following  4  models: 

(a)  apply  deterministic  boundary  conditions  obtained  from  the  steady  state  solution 

(b)  apply  stochastic  boundary  conditions  obtained  from  the  closed-form  steady  state  solution 

(c)  apply  approximate  stochastic  boundary  conditions,  ignoring  correlation  effects 

(d)  apply  approximate  stochastic  boundary  conditions,  accounting  for  the  correlation  between  the  boundaries 

Model  (a)  completely  ignores  boundary  condition  uncertainty.  Since  a  closed-form  solution  exists  for  this  par¬ 
ticular  problem  if  the  viscosity  is  modeled  by  a  single  random  variable,  we  can  get  a  fairly  accurate  modeling  of 
the  boundary  condition  uncertainty.  For  Model  (b),  we  applied  the  exact  PDF  obtained  for  the  boundary  conditions 
u(— a;max)  and  u(a;max)  to  the  random  held  model.  This  Model  (b)  still  represents  an  approximation  in  the  sense  that 
the  random  held  uncertainties  are  all  lumped  together  on  the  boundary. 

In  many  practical  cases  we  may  only  have  second  moment  information  about  the  boundary  condition  uncertainty 
(i.e.,  mean  value  and  variance).  It  then  makes  sense  to  add  a  random  perturbation  to  the  deterministic  boundary 
conditions  at  ±Xmax-  This  is  done  in  Models  (c)  and  (d).  Since  0  <  u{x)  <  1,  a  Gaussian  perturbation  at  the  boundary 
can  physically  not  be  justified.  Since  only  very  scant  information  may  be  available  regarding  the  uncertainty  at  the 
boundary  conditions,  we  used  the  following  basic  perturbation  model,  based  on  the  deterministic  boundary  conditions 
Eq.(10.2): 

I  U{-Xraa^)  =  0 -f  (5i 
\  u{Xms.^)  =  1  -  (52 

where  (5,  are  exponentially  distributed  random  variables,  with  variance  equal  to  the  variance  of  the  PDF  in  Model  (b). 
Model  (c)  ignores  the  correlation  effects  between  (5i  and  82',  while  these  are  taken  into  consideration  in  Model  (d).  The 
correlated  exponentially  distributed  random  variables  are  generated  using  the  Nataf  transformation  [8].  Figure  11.2 
compares  the  approximate  exponential  marginal  PDF  for  u(— a;max)  used  in  Models  (c)  and  (d)  with  the  PDF  used  in 
Model  (b). 

Figure  11.3  shows  that,  when  the  uncertainty  on  the  boundary  conditions  is  included,  the  difference  between 
the  deterministic  and  the  mean  stochastic  solution  is  three  times  larger  than  when  the  uncertainty  on  the  boundary 
conditions  is  not  included  (see  Figure  10.1).  For  practical  purposes,  this  difference  remains  small  though  (l^  =  1). 
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StDev[u(x)]  StDev[u(x)] 


Figure  a 


Figure  b 


Figure  c  Figure  d 

Fig.  11.1.  Comparison  o/StDev[t((a;)]  using  dijferent  boundary  condition  models  and  8193  grid  points. 


Note  that  the  mean  random  variable  solution  is  an  exact  result,  whereas  the  mean  random  held  solution  is  the  results 
of  10®  Monte  Carlo  simulations.  The  standard  error  on  the  Monte  Carlo  results  is  less  than  0.00006. 

Regarding  the  standard  deviation  of  u{x),  the  following  conclusions  can  be  drawn  from  the  analysis  (see  Figures 
11.1  and  11.4): 

•  Ignoring  the  uncertainty  of  the  boundary  conditions  grossly  underestimates  the  variability  in  the  interior 
domain  when  the  correlation  length  1^  is  substantially  shorter  than  the  integration  domain.  Note  that  this  will 
usually  be  the  case;  otherwise  there  is  no  reason  to  use  a  random  held  model. 

•  The  importance  of  an  accurate  estimation  of  the  correlation  between  random  variables  is  clear  from  a  com¬ 
parison  of  Models  (c)  and  (d).  If  the  correlation  length  is  long  and  substantial  correlation  is  present  between 
(5i  and  <52,  Models  (c)  and  (d)  give  drastically  different  results. 

•  Even  though  the  exponential  PDF  is  a  crude  approximation  of  the  boundary  uncertainty  (see  Figure  1 1.2),  the 
variance  of  Model  (d)  follows  the  trends  of  Model  (b)  remarkably  well  for  all  (^-values. 
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Fig.  11.2.  Comparison  between  approximate  exponential  and  exact  PDF  at  left-end  boundary,  computed  using  a  random  variable  model  for  p. 


-  Deterministic  -  Random  Field 

—  —  —  -  Random  Variable  -  Random  Field 


X 


Fig.  1 1.3.  Difference  between  deterministic,  mean  random  variable  and  mean  random  field  solution  (boundary  condition  Model  (b),  expo¬ 
nential  auto-correlation  with  1^  =  1). 


12.  Summary.  The  paper  describes  and  highlights  the  need  for  more  advanced  uncertainty  modeling.  The  dif¬ 
ferent  sources  of  uncertainty  and  their  treatment  are  reviewed.  The  steady-state  generalized  Burgers  equation  (GBE)  is 
solved  using  a  variety  of  stochastic  models.  First  a  random  variable  model  is  used.  The  differential  equation  is  solved 
deterministically  and  the  inherent  uncertainty  associated  with  the  viscosity  is  propagated  through  the  closed-form  so¬ 
lution.  The  exact  solutions  for  the  mean  and  standard  deviation  are  compared  with  Monte  Carlo  simulations  as  well  as 
with  first-  and  second-order  approximate  results.  The  second-order  results  are  found  to  be  significantly  more  accurate. 

Subsequently,  the  differential  equation  is  solved  stochastically  using  a  random  field  description  for  the  viscosity. 
The  viscosity  generally  depends  on  the  temperature  distribution  for  which  no  information  is  available  in  the  GBE. 
The  random  field  emulates  the  model  uncertainty  caused  by  these  missing  variables.  We  computed  results  using  both 
midpoint  and  locally  averaged  discretizations.  It  is  shown  how,  for  a  given  discretization  level,  both  methods  bracket 
the  true  variability.  They  provide  practical  upper  and  lower  bounds  to  the  variability.  The  impact  of  the  auto-correlation 
model  is  studied  as  well. 

The  use  of  a  random  field  model  for  the  viscosity  instead  of  a  single  random  variable  transforms  the  Burgers 
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Figure  a:  a;  =  —0.75  Figure  b:  a;  =  0 

Fig.  1 1 .4.  Impact  of  boundary  condition  uncertainty  on  StDev[t((a;)]/or  different  correlation  lengths  Ic- 


equation  into  a  stochastic  differential  equation.  The  presented  results  clearly  reveal  the  random  field  effect:  Var[u(0)] 
depends  on  both  the  correlation  length  1^  and  the  auto-correlation  model  R{^).  It  is  intuitively  clear  that  a  zero  variance 
at  a;  =  0,  as  obtained  from  the  random  variable  model,  is  unrealistic.  The  uncertainty  caused  by  the  missing  variable 
-  temperature  -  can  be  modeled  realistically  using  a  random  field  description  of  the  viscosity. 

In  the  last  section,  the  effects  of  boundary  condition  uncertainty  are  studied.  This  uncertainty  is  caused  by  the 
truncation  of  the  integration  domain  at  the  far  field  boundary  during  the  numerical  solution  phase  of  the  differential 
equations.  In  practice  only  limited  information  will  be  available  regarding  this  type  of  model  uncertainty.  Therefore 
we  compared  various  approximate  models.  It  can  be  concluded  that  ignoring  this  boundary  condition  uncertainty 
dramatically  underestimates  the  variance  of  the  velocity  u{x)  in  the  interior  of  the  domain. 
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